clear;
clc;
close all;
%% User input
fileExt = '.mat';
outputName = 'Figures';
dim = '3D'; %or '3D'
prop2Plot = 'extRad';
nStack = 5;
nConc = 3; %n of concentration
conc = {'_1_','_05','_025',};
len  = {'1K','5K'};
nColor = {'-r','-g','-b','-y','-k'};
nBins = 100; %for comparing different CDF
%% Load 
[file2Analyze,currentFolderName,outDir] = Load.Folder(fileExt,outputName);
data2Plot = Load.getData2Plot(file2Analyze,dim,prop2Plot);
fields = fieldnames(data2Plot);
nFile  = length(fields);
nCond = length(conc)*length(len);
%% GetCDF
for i = 1 : nFile
    
    currField = fields{i};
    allData = [];
    
    for j = 1:nStack
        currData  = data2Plot.(currField).(prop2Plot){j,:};
        [~,CCDF] = Plotting.getCDF(currData(:));
        data2Plot.(currField).CDF{j,1} = [CCDF.x CCDF.y];
        allData = [allData; currData(:)]; 
    end
         
    [~,CCDF] = Plotting.getCDF(allData);
    data2Plot.(currField).allCCDF{1,1} = [CCDF.x CCDF.y];
    data2Plot.(currField).allData{1,1} = allData;

end
% Here we combine S1 and S2 and we get the errors for error bar
for i = 1:nConc 
    idx    = contains(fields,conc{i});
    tmpField2Use = fields(idx);

    for j = 1:length(len)
        idx2 = contains(tmpField2Use,len{j});
        field2Use = tmpField2Use(idx2);
        
        S1  = data2Plot.(field2Use{1}).allData{1,1};
        S2  = data2Plot.(field2Use{2}).allData{1,1};
        combData = [S1(:);S2(:)];
        
        [~,CCDF] = Plotting.getCDF(combData);
        [midBin, Freq, LFreq] = Plotting.lnbin(CCDF.x,nBins);
        
        x = CCDF.x;
        
        meanSamples = zeros(nStack*2,length(midBin)-1);
        stdSamples = zeros(nStack*2,length(midBin)-1);
        minSamples = zeros(nStack*2,length(midBin)-1);
        maxSamples = zeros(nStack*2,length(midBin)-1);
        
        for k = 1:nStack
            dataS1 = data2Plot.(field2Use{1}).CDF{k,1};
            dataS2 = data2Plot.(field2Use{2}).CDF{k,1};
            y1 = dataS1(:,2);
            x1 = dataS1(:,1);
            y2 = dataS2(:,2);
            x2 = dataS2(:,1);
            for l = 1 : length(midBin)-1
                
                
                    meanSamples(2*k-1,l) = nanmean(y1(and(x1<=midBin(l+1),x1>midBin(l))));
                    min1S  = nanmin(y1(and(x1<=midBin(l+1),x1>midBin(l))));
                    max1S  = nanmax(y1(and(x1<=midBin(l+1),x1>midBin(l))));
                    if ~isempty (min1S)
                        minSamples(2*k-1,l) = min1S;
                        maxSamples(2*k-1,l) = max1S;
                    else
                        minSamples(2*k-1,l) = NaN;
                        maxSamples(2*k-1,l) = NaN;
                    end
            
                    meanSamples(2*k,l) = nanmean(y2(and(x2<=midBin(l+1),x2>midBin(l))));
                    min2S  = nanmin(y2(and(x2<=midBin(l+1),x2>midBin(l))));
                    max2S  = nanmax(y2(and(x2<=midBin(l+1),x2>midBin(l))));

                    if ~isempty (min2S)
                        minSamples(2*k,l) = min2S;
                        maxSamples(2*k,l) = max2S;
                    else
                        minSamples(2*k,l) = NaN;
                        maxSamples(2*k,l) = NaN;
                    end
             
            end
        end
        
        meanSamples = nanmean(meanSamples,1);
        minSamples = nanmin(minSamples,[],1);
        maxSamples = nanmax(maxSamples,[],1);
        %clean data
        midBin(isnan(meanSamples))=[];
        minSamples((isnan(meanSamples))) =[];
        maxSamples((isnan(meanSamples))) =[];
        meanSamples((isnan(meanSamples))) =[];
        midBin(end) = [];
        fieldName = ['cond',conc{i},'_',len{j}];
        combSample.(fieldName).midBin  = midBin(:);
        combSample.(fieldName).avg     = meanSamples(:);
        combSample.(fieldName).maxErr  = maxSamples(:)  - meanSamples(:);
        combSample.(fieldName).minErr  = meanSamples(:) - minSamples(:);
        combSample.(fieldName).allData = combData;
     
    end
end
%getCDF
%% Plotting Polymer Concentration
concLab = {'1mg/mL','0.5mg/mL','0.25mg/mL'};
%1K
idx    = contains(fields,'1K');
field2Use = fields(idx);
nfields = length(field2Use);

figure
%ALL curve on top of each other
for i = 1:nfields
    currField = field2Use{i};
    str(1:3) = {currField};
    test = cell2mat(cellfun(@contains,str,conc,'UniformOutput',0));
    currColor = nColor{test};
    for j = 1:nStack
        currCDF = data2Plot.(currField).CDF{j,1};
        subplot(1,3,1)
        hold on
        plot(currCDF(:,1),currCDF(:,2),currColor)
        ax = gca;
        ax.XScale = 'log';
        ax.YScale = 'log';
        axis square
        hold off
    end
end
names = fieldnames(combSample);
idx = contains(names,'1K');
field2Use = names(idx);
nfields = length(field2Use);

%S1 and S2 averaged
for i =1:nfields
    currField = field2Use{i};
    
    subplot(1,3,2)
    hold on
    
    err = [combSample.(currField).maxErr' ; combSample.(currField).minErr'];
    p = Plotting.shadedErrorBar(combSample.(currField).midBin,combSample.(currField).avg,err,'lineprops',nColor{i},'logScale',true);

    ax = gca;
    ax.XScale = 'log';
    ax.YScale = 'log';
    axis square
    hold off
    
    subplot(1,3,3)
    hold on
    h1 = histcounts(combSample.(currField).allData);
    plot(h1);
    axis square
    hold off
    
end
legend(concLab);
%5K
idx    = contains(fields,'5K');
field2Use = fields(idx);
nfields = length(field2Use);

figure
%ALL curve on top of each other
for i = 1:nfields
    currField = field2Use{i};
    str(1:3) = {currField};
    test = cell2mat(cellfun(@contains,str,conc,'UniformOutput',0));
    currColor = nColor{test};
    for j = 1:nStack
        currCDF = data2Plot.(currField).CDF{j,1};
        subplot(1,3,1)
        hold on
        plot(currCDF(:,1),currCDF(:,2),currColor)
        ax = gca;
        ax.XScale = 'log';
        ax.YScale = 'log';
        axis square
        hold off
    end
end
names = fieldnames(combSample);
idx = contains(names,'5K');
field2Use = names(idx);
nfields = length(field2Use);

%S1 and S2 averaged
for i =1:nfields
    currField = field2Use{i};
    
     subplot(1,3,2)
    hold on
    
    err = [combSample.(currField).maxErr' ; combSample.(currField).minErr'];
    p = Plotting.shadedErrorBar(combSample.(currField).midBin,combSample.(currField).avg,err,'lineprops',nColor{i},'logScale',true);

    ax = gca;
    ax.XScale = 'log';
    ax.YScale = 'log';
    axis square
    hold off
    
    subplot(1,3,3)
    hold on
     h1 = histcounts(combSample.(currField).allData);
    plot(h1);
    axis square
    hold off
end
legend(concLab);

%% Polymer Length dependency
conc = {'_1_','_05','_025',};
figure
for i = 1:nConc
    idx = contains(names,conc{i});
    field2Use = names(idx);
    for j = 1:length(field2Use)
        currField = field2Use{j};
        
        subplot(2,3,i)
        hold on
        err = [combSample.(currField).maxErr' ; combSample.(currField).minErr'];
        p = Plotting.shadedErrorBar(combSample.(currField).midBin,combSample.(currField).avg,err,'lineprops',nColor{j},'logScale',true);

        ax = gca;
        ax.XScale = 'log';
        ax.YScale = 'log';
        axis square
        hold off
        title(concLab{i});
        subplot(2,3,nConc+i)
        hold on
         h1 = histcounts(combSample.(currField).allData);
    plot(h1);
        axis square
        hold off
    end
end